/////// 求解二维Possion或者Laplace方程,  形式为Delat(P) = F, 右端项F可以等于零 

#include <time.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <hip/hip_runtime.h>

#include "kernel.h"
#include "util.h"

#include "main.h"

#define REAL double

#define USE_CUDA

int nx, ny, n_level, p_relax;
REAL pi, h, dt, xleft, xright, yleft, yright, Re, We, shear_rate, velocity_boundary;
REAL 	**bd_unit_normal, *bd_curvature;  // area_correct() use these variable
REAL epsilon;

//////////////  测试核函数！！！！！！！  
int main(int argc, char** argv)
{
	char str1[8], str2[8], str3[8];
	sprintf(str1, "./data");
	sprintf(str2, "./uvp");
	sprintf(str3, "./bdpo");
	if (!access(str1, 0) == 0) mkdir(str1, 0777);
	if (!access(str2, 0) == 0) mkdir(str2, 0777);
	if (!access(str3, 0) == 0) mkdir(str3, 0777);


	extern int nx, ny, n_level, p_relax;
	extern REAL pi, h, dt, xleft, xright, yleft, yright, Re, We, shear_rate, velocity_boundary, epsilon;
	extern REAL **bd_unit_normal, *bd_curvature;

	FILE * fp, *my;
	char buf[64];

	clock_t start, end;
	REAL elapsed;

	pi = 4.0*atan(1.0);  // 圆周率
	p_relax = 20;

	////////////////////////////
	int nx = gnx;
	int ny = gny;
		
	xleft  = 0.0;
	yleft  = 0.0;

	yright = 8.0;
	xright = (REAL)(nx / ny)*yright;
	
	REAL h = (xright - xleft) / (REAL)nx;
	REAL h2 = h*h;

	REAL endT = 1.0;
	REAL dt = 0.5* h2;
	
	int max_it = (int)(endT / dt + 0.001);
	int ns = (int)(max_it / 5 + 0.001);

	REAL radi = 0.5; // 液滴半径
	int nt = (int)(2 * pi * radi / (h / 5.0) + 0.001);

	REAL shear_rate = 1.0; //velocity                // 初始上下边界剪切速度
	REAL Ca = 0.3;                          // 流体力学中的Capillary number
	REAL Re = 20.0;                          // 流体力学中的雷诺数

	REAL kappa = shear_rate / (Ca*Re); //mu=1          // 流体力学中表面张力系数

	//REAL kappa = 0.005;

	/////////////////////////////////////////////////

	n_level = (int)(log(ny) / log(2) - 0.9);
	if (n_level > (int)(log(nx) / log(2) - 0.9))
		n_level = (int)(log(nx) / log(2) - 0.9);

	printf("n_level = %d\n", n_level);

	char bufferremarks[128];
	sprintf(bufferremarks, "remarks.m");
	my = fopen(bufferremarks, "w");
	fprintf(my, "nx          = %d\n", nx);
	fprintf(my, "ny          = %d\n", ny);
	fprintf(my, "xright      = %f\n", xright);
	fprintf(my, "yright      = %f\n", yright);

	fprintf(my, "max_it      = %d\n", max_it);
	fprintf(my, "dt          = %f\n", dt);
	fprintf(my, "nt          = %d\n", nt);

	fprintf(my, "out_count   = %d\n", (int)(max_it / ns + 0.001));
	fprintf(my, "h           = %f\n", h);
	fprintf(my, "Re          = %f\n", Re);
	fprintf(my, "ra          = %f\n", radi);
	fprintf(my, "total_count = %d\n", max_it / ns);

	fprintf(my, "kappa       = %f\n", kappa);
	fprintf(my, "Ca          = %f\n", Ca);
	fprintf(my, "n_level     = %d\n\n", n_level);
	fclose(my);

	//fclose(my);
	//printf("123456789987654321\n");
	
	REAL **bdpo = dmatrix(0, nt + 1, 1, 2);  // 界面边界点坐标数组
	REAL *ds    = dvector(0, nt + 1);        // 边界点相邻两点之间的距离

	REAL **fxx = dmatrix(0, nx, 1, ny);     // 整体背景网格下的流体力, x方向
	REAL **fyy = dmatrix(1, nx, 0, ny);     // 整体背景网格下的流体力, y方向

	REAL **pmat = dmatrix(1, nx, 1, ny); 
	REAL **fmat = dmatrix(1, nx, 1, ny);

	initial(pmat, fmat, 1, nx, 1, ny, xright, yright, bdpo, nt, radi);  // 初始化二维矩阵p和f
	//printf("123456789987654321\n");
    // 由界面边界点信息,得到界面上的边界力,bdforcex和bdforcey,用delta函数转化到整个计算域的流体力fxx和fyy	
	ibm_force(bdpo, ds, fxx, fyy, h,kappa, nx, ny, nt);
	// 对fxx和fyy求散度,即我们讲过的用delta_dot作用于fxx和fyy. 输出就是和pmat维数相同的二维矩阵了 
	div_uv(fxx, fyy, fmat, nx, ny);

	// 到这里,我们有了维数相同的二维矩阵pmat和fmat, 然后转化为一维数组用GPU计算,这个Possion方程....
	// 下面的是我们已经实现的并行代码, 二维数组转一维等初始化在函数cpu_relax_p()中， 然后调用函数relax_p_kernelLauncher(),
	// 包装函数relax_p_kernelLauncher()调用核函数relax_p_kernel()在GPU中实际计算!!!
	//gpu_relax_p(pmat, fmat, nx, ny, nx);

	MG_Poisson(pmat, fmat, nx, ny);

	// 输出二维矩阵结果到文件中
	sprintf(buf, "pmat.m");
	fp = fopen(buf, "w");
	for (int i = 1; i <= nx; i++)
	{
		for (int j = 1; j <= ny; j++)
		{
			fprintf(fp, "  %5.6f", pmat[i][j]);
		}
		fprintf(fp, "\n");
	}
	fclose(fp);
	//printf("123456789987654321\n");
	sprintf(buf, "fmat.m");
	fp = fopen(buf, "w");
	for (int i = 1; i <= nx; i++)
	{
		for (int j = 1; j <= ny; j++)
		{
			fprintf(fp, "  %5.6f", fmat[i][j]);
		}
		fprintf(fp, "\n");
	}
	fclose(fp);


	// 或者写函数调用函数print_mat()    
	// 最后在matlab中画图

	free_dmatrix(pmat, 1, nx, 1, ny);
	free_dmatrix(fmat, 1, nx, 1, ny);
	system("pause");

	return 0;
}


void MG_Poisson(REAL **p, REAL **f, int nx, int ny)     // 多重网格方法求解压力泊松方程
{
	int    max_it = 4000, it_mg = 1;
	REAL tol = 1.0e-6, resid = 10.0;

	REAL **sor = dmatrix(1, nx, 1, ny);
	REAL **workv = dmatrix(1, nx, 1, ny);

	mat_copy(workv, p, 1, nx, 1, ny);

	
	while (it_mg <= max_it && resid >= tol)
	{

		vcycle_uv(p, f, nx, ny, nx, 1);

		pressure_update(p, nx, ny);


		ijloop
		{
			sor[i][j] = workv[i][j] - p[i][j];
		}

		resid = mat_max(sor, 1, nx, 1, ny);
		mat_copy(workv, p, 1, nx, 1, ny);

		it_mg++;
	}


	printf("Mac iteration = %d  residual = %16.15f \n", it_mg, resid);

	free_dmatrix(sor,   1, nx, 1, ny);
	free_dmatrix(workv, 1, nx, 1, ny);

	return;
}

void vcycle_uv(REAL **uf, REAL **ff, int nxf, int nyf, int nx, int ilevel)
{
	extern int n_level;

	gpu_relax_p(uf, ff, nxf, nyf);
	//cpu_relax_p(uf, ff, nxf, nyf);

	//printf("MG_debug\n");
	//system("pause");
	//exit(99);

	//if (ilevel < n_level)
	//{
	//	int nxc, nyc;
	//	REAL **rf, **uc, **fc;

	//	nxc = nxf / 2;
	//	nyc = nyf / 2;

	//	rf = dmatrix(1, nxf, 1, nyf);
	//	uc = dmatrix(1, nxc, 1, nyc);
	//	fc = dmatrix(1, nxc, 1, nyc);

	//	residual_p(rf, uf, ff, nxf, nyf);

	//	restrict1(rf, fc, nxc, nyc);
	//	zero_matrix(uc, 1, nxc, 1, nyc);
	//	vcycle_uv(uc, fc, nxc, nyc, nx, ilevel + 1);

	//	prolong(uc, rf, nxc, nyc);
	//	mat_add(uf, uf, rf, 1, nxf, 1, nyf);

	//	//gpu_relax_p(uf, ff, nxf, nyf);
	//	cpu_relax_p(uf, ff, nxf, nyf);

	//	free_dmatrix(rf, 1, nxf, 1, nyf);
	//	free_dmatrix(uc, 1, nxc, 1, nyc);
	//	free_dmatrix(fc, 1, nxc, 1, nyc);
	//}

}


void restrict1(REAL **u_fine, REAL **u_coarse, int nxt, int nyt)
{
	ijloopt
	{
		u_coarse[i][j] = 0.25*(u_fine[2 * i - 1][2 * j - 1]
		+ u_fine[2 * i - 1][2 * j]
			+ u_fine[2 * i][2 * j - 1]
			+ u_fine[2 * i][2 * j]);
	}
}


void prolong(REAL **u_coarse, REAL **u_fine, int nxt, int nyt)
{
	ijloopt
	{
		u_fine[2 * i - 1][2 * j - 1] = u_fine[2 * i - 1][2 * j] =
		u_fine[2 * i][2 * j - 1] = u_fine[2 * i][2 * j] = u_coarse[i][j];
	}
}

void residual_p(REAL **r, REAL **u, REAL **f, int nxt, int nyt)
{
	laplace_p(u, r, nxt, nyt);
	mat_sub(r, f, r, 1, nxt, 1, nyt);
}


void laplace_p(REAL **p, REAL **lap_p, int nxt, int nyt)
{
	REAL **dpdx, **dpdy;

	dpdx = dmatrix(0, nxt, 1, nyt);
	dpdy = dmatrix(1, nxt, 0, nyt);

	grad_p(p, dpdx, dpdy, nxt, nyt);
	div_uv(dpdx, dpdy, lap_p, nxt, nyt);

	free_dmatrix(dpdx, 0, nxt, 1, nyt);
	free_dmatrix(dpdy, 1, nxt, 0, nyt);
}


void grad_p(REAL **p, REAL **dpdx, REAL **dpdy, int nxt, int nyt)
{
	extern REAL xright;

	REAL ht;

	ht = xright / (REAL)nxt;

	i0jloopt
	{
		if (i == 0)
		dpdx[0][j] = (p[1][j] - p[nxt][j]) / ht;

		else if (i == nxt)
			dpdx[nxt][j] = (p[1][j] - p[nxt][j]) / ht;

		else
			dpdx[i][j] = (p[i + 1][j] - p[i][j]) / ht;
	}

		// if(nxt==gnx)
	{
		ij0loopt
		{
			if (j == 0)
			dpdy[i][0] = 0.0;
			else if (j == nyt)
				dpdy[i][nyt] = 0.0;
			else
				dpdy[i][j] = (p[i][j + 1] - p[i][j]) / ht;
		}
	}
	// else
	// {  
	//    ij0loopt 
	//    {	  
	//  	  if (j==0)
	//            dpdy[i][0] = (p[i][1] - p[i][nyt])/ht; 
	//  	  else if (j==nyt)
	//  		  dpdy[i][nyt] = (p[i][1] - p[i][j])/ht; 
	//        else
	//            dpdy[i][j] = (p[i][j+1] - p[i][j])/ht;
	//     }  
	//  }
}



void pressure_update(REAL** a, int nx, int ny)
{
	int i, j;
	REAL ave = 0.0;

	for (i = 1; i <= nx; i++)
		for (j = 1; j <= ny; j++)
		{
			ave = ave + a[i][j];
		}

	ave /= (nx + 0.0) * (ny + 0.0);
	for (i = 1; i <= nx; i++)
		for (j = 1; j <= ny; j++)
		{
			a[i][j] -= ave;
		}
	return;
}


void  gpu_relax_p(REAL **pmatf, REAL **fmatf, int nxf, int nyf)
{
	///////////////////////////////////////////////////////////////////////////

	extern int p_relax;

	extern REAL xright;

	REAL ht = xright / (REAL)nxf;

	//REAL** ap = dmatrix(0, nxf + 1, 0, nyf + 1);	
	//REAL** af = dmatrix(0, nxf + 1, 0, nyf + 1);
	//zero_matrix(ap, 0, nxf + 1, 0, nyf + 1);
	//zero_matrix(af, 0, nxf + 1, 0, nyf + 1);

	//mat_copy(ap, pmatf, 1, nxf, 1, nyf);
	//mat_copy(af, fmatf, 1, nxf, 1, nyf);

	//augmen_p(ap, pmatf, nxf, nyf);
	//augmen_p(af, fmatf, nxf, nyf);

	// p和f开辟一维内存数组
	REAL  *pvec = dvector(0, (nxf)* (nyf) - 1);
	REAL  *fvec = dvector(0, (nxf)* (nyf) - 1);

	// 开辟对应的一维显存
	REAL* d_pvec, *d_fvec;
	hipMalloc(&d_pvec, (nxf) * (nyf) * sizeof(REAL));
	hipMalloc(&d_fvec, (nxf) * (nyf) * sizeof(REAL));

	mat2vec(pvec, pmatf, 1, nxf, 1, nyf);       // 二维矩阵转化为对应的一维向量

	mat2vec(fvec, fmatf, 1, nxf, 1, nyf);      // 二维矩阵转化为对应的一维向量

	// 对应一维向量拷贝到显存
	hipMemcpy(d_pvec, pvec, (nxf) * (nyf) * sizeof(REAL), hipMemcpyHostToDevice);
	hipMemcpy(d_fvec, fvec, (nxf) * (nyf) * sizeof(REAL), hipMemcpyHostToDevice);

	// 调用GPU核心计算
	for (int i = 1; i <= p_relax; i++)
	{
		//printf("iter = %d\n", i);
		///// test kernel ////
		relax_p_kernelLauncher(d_pvec, d_fvec, nxf, nyf, ht);
	}

	//// 对应一维向量拷贝回主机内存
	hipMemcpy(pvec, d_pvec, (nxf) * (nyf) * sizeof(REAL), hipMemcpyDeviceToHost);
	hipMemcpy(fvec, d_fvec, (nxf) * (nyf) * sizeof(REAL), hipMemcpyDeviceToHost);

	hipFree(d_pvec);
	hipFree(d_fvec);

	//对应一维向量转化为二维矩阵
	vec2mat(pmatf, pvec, 1, nxf, 1, nyf);
	vec2mat(fmatf, fvec, 1, nxf, 1, nyf);

	//mat_copy(pmatf, ap, 1, nxf, 1, nyf);
	//mat_copy(fmatf, af, 1, nxf, 1, nyf);

	free_dvector(pvec, 0, (nxf)*(nyf) - 1);
	free_dvector(fvec, 0, (nxf)*(nyf) - 1);

	//free_dmatrix(ap, 0, nxf + 1, 0, nyf + 1);
	//free_dmatrix(af, 0, nxf + 1, 0, nyf + 1);
}


void initial(REAL **pmat, REAL **fmat, int nxl, int nxr, int nyl, int nyr, REAL xright, REAL yright, REAL **bdpo, int nt, REAL radi)
{
	extern REAL pi;

	for (int i = nxl; i <= nxr; i++)
		for (int j = nyl; j <= nyr; j++)
		{
			pmat[i][j] = 0.0;     // 初始化        
		}

	for (int i = nxl; i <= nxr; i++)
		for (int j = nyl; j <= nyr; j++)
		{
			fmat[i][j] = 0.0;
		}

	for (int i = 1; i <= nt; i++)
	{
		bdpo[i][1] = (radi*cos(2.0*pi*i / nt) + 0.5*(0.0 + xright));     // 产生初始的液滴边界，圆形，ra为半径，ibm方法中，phi是圆周的边界点坐标
		bdpo[i][2] = (radi*sin(2.0*pi*i / nt) + 0.5*(0.0 + yright));
	}

	bdpo[0][1] = bdpo[nt][1];                // 因为是数组表示，保证首尾相连
	bdpo[nt + 1][1] = bdpo[1][1];
	bdpo[0][2] = bdpo[nt][2];
	bdpo[nt + 1][2] = bdpo[1][2];

}

//////////////////////////////////////////////////////////
void augmen_p(REAL** ap, REAL** p, int nxt, int nyt)
{
	int i, j;

	mat_copy(ap, p, 1, nxt, 1, nyt);

	for (i = 0; i <= nxt; i++)
	{
		//ap[i][0] = ap[i][1];
		//ap[i][nyt + 1] = ap[i][nyt];

		//测试
		ap[i][0] = 0.0;
		ap[i][nyt + 1] = 0.0;
	}

	for (j = 0; j <= nyt; j++)
	{
		//ap[0][j] = ap[nxt][j];
		//ap[nxt + 1][j] = ap[1][j];

		//测试
		ap[0][j] = 0.0;
		ap[nxt + 1][j] = 0.0;
	}
}

 void cpu_relax_p(REAL **p, REAL **f, int nxt, int nyt)
 {
     extern int p_relax;
     extern REAL xright;
 
     int    iter;
     REAL ht, ht2, coef, src;
 
     ht = xright / (REAL) nxt;
     ht2 = pow(ht, 2);
 

     REAL** ap = dmatrix( 0, nxt + 1, 0, nyt + 1);
 
     zero_matrix(ap, 0, nxt + 1, 0, nyt + 1);

     for (iter = 1; iter <= p_relax; iter++) 
     {      
         augmen_p(ap, p, nxt, nyt);
 
         ijloopt 
         {
 		    coef = -4.0/ht2;
 		    //src = f[i][j];
             src = -(ap[i + 1][j] + ap[i - 1][j] + ap[i][j + 1] + ap[i][j - 1]) / ht2 + f[i][j];

             p[i][j] = src / coef;    

			/* if (i > 9 && i < 12)
			 {
				 printf("idx=%d  f[%d][%d]=%lf\n",i*32+j,i,j,f[i][j]);
			 }*/
 	    }     
     }
 
     free_dmatrix(ap, 0, nxt + 1, 0, nyt + 1);
 }
///////////////////////////////////////////////////////////////////////////////////////////////


//void cpu_relax_p(REAL** p, REAL** f, int nxt, int nyt)
//{
//    extern int    p_relax;
//    extern REAL xright;
//
//    int i, j, iter;
//    REAL h2, coef, src;
//
//    h2 = pow(xright / (REAL)nxt, 2);
//
//    for (iter = 1; iter <= p_relax; iter++) 
//    {
//        ijloopt
//        {
//            src = f[i][j];
//
//            if (i > 1 && i < nxt) 
//            {
//                src -= (p[i + 1][j] + p[i - 1][j]) / h2;
//                coef = 2.0 / h2;
//            }
//            else if (i == 1) 
//            {
//                src -= (p[i + 1][j] + p[nxt][j]) / h2;
//                coef = 2.0 / h2;
//            }
//            else 
//            {
//                src -= (p[1][j] + p[i - 1][j]) / h2;
//                coef = 2.0 / h2;
//            }
//
//            if (j > 1 && j < nyt) 
//            {
//                src -= (p[i][j + 1] + p[i][j - 1]) / h2;
//                coef += 2.0 / h2;
//            }
//            else if (j == 1) 
//            {
//                src -= p[i][2] / h2;
//                coef += 1.0 / h2;
//            }
//            else 
//            {
//                src -= p[i][nyt - 1] / h2;
//                coef += 1.0 / h2;
//            }
//
//            p[i][j] = -src / coef; 
//        }
//    }
//}



//void cpu_relax_p(REAL** p, REAL** f, int nxt, int nyt)
//{
//   extern int p_relax;
//   extern REAL xright;
//
//   int i, j, iter;
//   REAL ht, ht2, coef, src;
//
//   ht = xright / (REAL)nxt;
//   ht2 = pow(ht, 2);
//
//   for (iter = 1; iter <= p_relax; iter++)
//   {
//       ijloopt
//       {
//          coef = -2.0 / ht2;
//          src = f[i][j];
//
//          if (i == 1)
//          {
//             src -= (p[2][j] + p[nxt][j]) / ht2;
//          }
//          else if (i == nxt)
//          {
//             src -= (p[1][j] + p[nxt - 1][j]) / ht2;
//          }
//          else
//          {
//             src -= (p[i + 1][j] + p[i - 1][j]) / ht2;
//          }
//
//          if (j == 1)
//          {
//             src -= (p[i][2]) / ht2;
//             coef += -1.0 / ht2;
//          }
//          else if (j == nyt)
//          {
//             src -= (p[i][nyt - 1]) / ht2;
//             coef += -1.0 / ht2;
//          }
//          else
//          {
//             src -= (p[i][j + 1] + p[i][j - 1]) / ht2;
//             coef += -2.0 / ht2;
//           }
//
//           p[i][j] = src / coef;
//       }
//   }
//}

/////////////////////////////////////////////////////////////////////////////



//////////////////////////////////ibm force//////////////////////////////////

REAL delta(REAL x)
{
	REAL value;

	if (fabs(x) > 2.0)
		value = 0.0;
	else if (fabs(x) < 1.0)
		value = (3.0 - 2.0*fabs(x) + sqrt(1.0 + 4.0*fabs(x) - 4.0*x*x)) / 8.0;
	else
		value = 0.5 - (3.0 - 2.0*(2.0 - fabs(x)) + sqrt(1.0 + 4.0*(2.0 - fabs(x)) - 4.0*pow(2.0 - fabs(x), 2))) / 8.0;

	return value;
}

// ibm方法中计算边界力， 转化流场体积力，在这个函数中是fxx和fyy，
// 注意fxx和fyy应该是全局变量，这两个方向的力fxx 和fyy 最后要加到NS返程的右端作为源项，模拟边界力作用
// nt是表示界面的边界点的个数
void ibm_force(REAL** bd_posi, REAL* ds, REAL **fxx, REAL **fyy, REAL h, REAL kappa, int nx, int ny, int nt) 
{
	int i, j, k, ia, ja;
	REAL d, d1,  bd_forcex, bd_forcey;

	REAL h2 = h*h;

	REAL **fx = dmatrix(0, nx + 1, 0, ny + 1);
	REAL **fy = dmatrix(0, nx + 1, 0, ny + 1);


	zero_matrix(fx, 0, nx + 1, 0, ny + 1);
	zero_matrix(fy, 0, nx + 1, 0, ny + 1);

	for (k = 1; k <= nt + 1; k++) {
		ds[k] = sqrt(pow(bd_posi[k][1] - bd_posi[k - 1][1], 2) + pow(bd_posi[k][2] - bd_posi[k - 1][2], 2));
	}
	ds[0] = ds[nt];

	for (k = 1; k <= nt; k++) 
	{
		ia = (int)(bd_posi[k][1] / h + 0.499999) - 1;
		ja = (int)(bd_posi[k][2] / h + 0.499999) - 1;

		bd_forcex = kappa * ((bd_posi[k + 1][1] - bd_posi[k][1]) / ds[k + 1] - (bd_posi[k][1] - bd_posi[k - 1][1]) / ds[k]);

		bd_forcey = kappa * ((bd_posi[k + 1][2] - bd_posi[k][2]) / ds[k + 1] - (bd_posi[k][2] - bd_posi[k - 1][2]) / ds[k]);


		for (i = ia; i <= ia + 3; i++) 
		{
			d1 = delta((REAL)i - 0.5 - bd_posi[k][1] / h) / h2;
			for (j = ja; j <= ja + 3; j++) 
			{
				d = d1 * delta((REAL)j - 0.5 - bd_posi[k][2] / h);
				fx[i][j] += d * bd_forcex;
				fy[i][j] += d * bd_forcey;
			}
		}
	}

	ijloop
	{
		if (i == nx)
	       fxx[0][j] = fxx[nx][j] = 0.5 * (fx[nx][j] + fx[1][j]);
		else
		   fxx[i][j] = 0.5 * (fx[i][j] + fx[i + 1][j]);

	    if (j == ny)
		   fyy[i][0] = fyy[i][ny] = 0.5 * (fy[i][ny] + fy[i][1]);
	    else
		   fyy[i][j] = 0.5 * (fy[i][j] + fy[i][j + 1]);
	}

	free_dmatrix(fx, 0, nx + 1, 0, ny + 1);
	free_dmatrix(fy, 0, nx + 1, 0, ny + 1);
}

///////////////////////////////////////////////////////////////////////////////////////////

void div_uv(REAL **tu, REAL **tv, REAL **divuv, int nxt, int nyt)
{
	extern REAL xright;

	REAL ht = xright / (REAL)nxt;

	ijloopt
	{
		divuv[i][j] = (tu[i][j] - tu[i - 1][j] + tv[i][j] - tv[i][j - 1]) / ht;
	}
}



//////////////////////////////////////////////////////////////////////////////////////////////////////

//int main(int argc, char** argv)
//{
//
//	extern int n_level;
//
//	FILE * fp, *my;
//	char buf[64];
//
//	clock_t start, end;
//	REAL elapsed;
//
//	BC bc;
//
//	int nx = gnx;
//	int ny = gny;
//
//	n_level = (int)(log(ny) / log(2) - 0.9);
//	if (n_level > (int)(log(nx) / log(2) - 0.9))
//		n_level = (int)(log(nx) / log(2) - 0.9);
//
//
//#ifdef USE_CUDA
//	printf("CUDA computing....\n");
//#else
//	printf("CPU computing....\n");
//#endif
//
//	char bufferremarks[128];
//	sprintf(bufferremarks, "remarks.m");
//	my = fopen(bufferremarks, "w");
//	fprintf(my, "nx = %d\n", nx);
//	fprintf(my, "ny = %d\n", ny);
//	fclose(my);
//
//	REAL **pmat = dmatrix(1, nx, 1, ny);
//	REAL **fmat = dmatrix(1, nx, 1, ny);
//
//	REAL *pvec = dvector(0, nx * ny - 1);
//	REAL *fvec = dvector(0, nx * ny - 1);
//
//	start = clock();
//
//	initial(pmat, pmat, 1, nx, 1, ny);  // 初始化输入数据
//
//	mat2vec(pvec, pmat, 1, nx, 1, ny);              // 二维矩阵转化为对应的一维向量
//	mat2vec(fvec, fmat, 1, nx, 1, ny);              // 二维矩阵转化为对应的一维向量
//
//
//#ifdef USE_CUDA
//													// 开辟对应的一维显存
//	REAL* d_pvec, *d_fvec;
//	hipMalloc(&d_pvec, nx * ny * sizeof(REAL));
//	hipMalloc(&d_fvec, nx * ny * sizeof(REAL));
//
//	// 对应一维向量拷贝到显存
//	hipMemcpy(d_pvec, pvec, (nx * ny) * sizeof(REAL), hipMemcpyHostToDevice);
//	hipMemcpy(d_fvec, fvec, (nx * ny) * sizeof(REAL), hipMemcpyHostToDevice);
//#endif // USE_CUDA
//
//	// 调用GPU核心计算
//	for (int i = 1; i <= ITERS_PER_RENDER; i++)
//	{
//		printf("iter = %d\n", i);
//
//#ifdef USE_CUDA
//		kernelLauncher(d_pvec, d_fvec, nx, ny, nx);
//#else
//		cpu_possion(pmat, fmat, nx, ny);
//#endif
//	}
//
//#ifdef USE_CUDA
//	// 对应一维向量拷贝回主机内存
//	hipMemcpy(pvec, d_pvec, nx * ny * sizeof(REAL), hipMemcpyDeviceToHost);
//	hipMemcpy(fvec, d_fvec, nx * ny * sizeof(REAL), hipMemcpyDeviceToHost);
//
//	hipFree(d_pvec);
//	hipFree(d_fvec);
//
//	//对应一维向量转化为二维矩阵
//	vec2mat(pmat, pvec, 1, nx, 1, ny);
//	vec2mat(fmat, fvec, 1, nx, 1, ny);
//#endif
//
//	end = clock();
//	elapsed = ((REAL)(end - start)) / CLOCKS_PER_SEC;
//	printf("Time elapsed %f\n", elapsed);
//
//
//
//	// 输出二维矩阵结果到文件中
//	sprintf(buf, "pmat.m");
//	fp = fopen(buf, "w");
//	for (int i = 1; i <= nx; i++)
//	{
//		for (int j = 1; j <= ny; j++)
//		{
//			fprintf(fp, "  %5.6f", pmat[i][j]);
//		}
//		fprintf(fp, "\n");
//	}
//	fclose(fp);
//
//	// 或者写函数调用函数print_mat()    
//	// 最后在matlab中画图
//
//	free_dmatrix(pmat, 1, nx, 1, ny);
//	free_dmatrix(fmat, 1, nx, 1, ny);
//
//	free_dvector(pvec, 0, nx*ny - 1);
//	free_dvector(fvec, 0, nx*ny - 1);
//
//	system("pause");
//
//	return 0;
//}

////////////////////////////////////////////////////////////////////////////////////////



// useful test functions!!!!
//REAL gpu_res = 0;
//dotLauncher(&gpu_res, pvec, fvec, nx*ny);
//printf("gpu result = %d\n", gpu_res);

//sprintf(buf, "pvec.m");
//fp = fopen(buf, "w");
//for (int i = 0; i < nx2 * ny2; i++)
//{
//	fprintf(fp, "  %5.6f", pvec[i]);

//	fprintf(fp, "\n");
//}
//fclose(fp);
//exit(0);

//// 输出二维矩阵结果到文件中
//sprintf(buf, "pmat.m");
//fp = fopen(buf, "w");
//for (int i = 1; i <= nx; i++)
//{
//	for (int j = 1; j <= ny; j++)
//	{
//		fprintf(fp, "  %5.6f", pmat[i][j]);
//	}
//	fprintf(fp, "\n");
//}
//fclose(fp);
